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3 ■ ABSTRACT 
(N- 

! The Haumea family is currently the only identified collisional family in the 

Ph " Kuiper belt. We numerically simulate the long-term dynamical evolution of the 

family to estimate a lower limit of the family's age and to assess how the popula- 
tion of the family and its dynamical clustering are preserved over Gyr timescales. 
^ I We find that the family is not younger than 100 Myr, and its age is at least 1 

^ \ Gyr with 95% confidence. We find that for initial velocity dispersions of 50 — 400 

^ I ms~^, approximately 20 — 45% of the family members are lost to close encounters 

\ with Neptune after 3.5 Gyr of orbital evolution. We apply these loss rates to two 

CN I proposed models for the formation of the Haumea family, a graze-and-merge type 

Q>^ I collision between two similarly sized, differentiated KBOs or the collisional dis- 

^ I ruption of a satellite orbiting Haumea. For the graze-and-merge collision model, 

I we calculate that > 85% of the expected mass in surviving family members within 

\0 \ 150 ms~^ of the collision has been identified, but that one to two times the mass 

^ I of the known family members remains to be identified at larger velocities. For the 

^ I satellite-break-up model, we estimate that the currently identified family mem- 

I bers account for ~ 50% of the expected mass of the family. Taking observational 

^ ■ incompleteness into account, the observed number of Haumea family members is 

I consistent with either formation scenario at the la level, however both models 

predict more objects at larger relative velocities (> 150 ms~^) than have been 
identified. 



Introduction 



The Haumea (2003 ELgi) collisional family was discovered by lBrown et all (120071 ) who 
noted that Haumea and five other Kuiper belt objects (KBOs) shared a spectral feature that 
is indicative of nearly pure water ice on the sur faces of the bodies. These s ix KBOs, along 
with four additional family members identified by lSchaller and Brown! fj2008l ). ISnodgrass et al. 
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( 120101 ). and iRagozzine and Brown! (120071 ). can all be dynamically linked to Haumea, and 
there do not appear to be any dynamically unrelated KBOs that share this spectral fea- 
ture. Aside from being spectrally linked to these other KBOs, Haumea itself shows signs 
of its coUisiona l past. Despite having a nearly pure water ice surface, Haumea's density is 
~ 2.6 g cm~^ ( iRabinowitz et a/.l 120061). which i s higher than expected for typical assumed 
ice/rock ratios in the Kuiper belt ( iBrownl 120081 ): one way to achieve this higher density is 
to have a catastrophic collision between a differentiated proto-Haumea and a nother KBQ 



in wh ich proto-Haumea loses a substantial fraction of its water ice mantle (IBrown et al. 
20071). This scenario is supported by the prese nce of at least two water ice satellites 



( iBarkume et al\ l2006l : IRagozzine and Brown! |2009| ) . Haumea also has a n elongated shape 
and a very short spin period of ~ 4 hours that is unlikely to be primordial ( IRabinowitz et al 
20061 : Ihacerda and Jewit"tll2007h . 



Ragozzine and Brown! (12007! ) examined the dynamical connections between the identi- 



fied Haumea family members. These connections are made by first estimating the orbit of 
the center of mass of the colliding bodies, and then estimating the ejection velocities of each 
family member relative to the collision's center of mass. The ejection velocity is given by 

Av = v-Vr^, d) 



where Vcm is the estimated collision's center-of-mass velocity. Because Haumea is by far 
the largest remnant from the collision, its orbit immediately after the collision should 
have nearly coincided with the center-of-mass orbit. However, Haumea is currently lo- 
cated at the boundary of the 12:7 mean motion resonance (MMR) with Neptune; over long 
timescales, the chaotic zone of this resonance causes a random walk of the proper elements 
such that Haumea's current orbit may be significantly distant from its post-collision orbit. 
Ragozzine and Brownl (!2007l ) estimate the center-of-mass collision orbit by minimizing the 
sum of the relative speeds of all family members, assuming that Haumea's semimajor axis 
and its Tisserand parameter with respect to Neptune are both conserved during its chaotic 
evolution; they then use Haumea's present distance from the collision's center-of-mass orbit, 
together with a calculation of its chaotic diffusion rate, to estimate the age of the collisional 
family to be 3.5 ± 2 Gy. Given the exceedingly low collision probabilities for objects large 
enough to form the Haumea family in the current Kuiper belt, the family is likely to be old. 
However, the family probably cannot have formed in the primordial, much more massive 
Kuiper belt, because whatever caused the mass of the Kuiper belt to be depleted (by an 
estimated 2 or 3 orders of rnagnit ude) would have also destroyed the dynamical coherence 
of the family (ILevison et a/.!!2008l ). The high inclination (~ 27°) of the family also argues 
against a primordial origin, because such large inclinations are probably products of the 
excitation and mass depletion of the Kuiper belt. Thus, it appears that the Haumea family- 
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forming collision occurred near the end of the primordial, high-mass phase of the Kuiper 
beh. 



Several of the largest KBOs show evidence of their coUisional past (see review by iBrown 



( I2OO8I )). but the Haumea family is the only collisional family that has been identified in the 
Kuiper belt. The dynamical connections between the members of the family allow us to 
place some constraints on the type of collision that formed the family and also constrain the 
age of the family as being old, but probably not primordial. These characteristics make the 
Haumea family an excellent probe of the collisional environment in the Kuiper belt following 
the excitation and mass depletion event; understanding the type of collision that created the 
family (especially the relative sizes and speeds of the impactor and target) would provide 
valuable insight into the size and orbital distribution of th e Kuiper belt at the time of the 
colhsion (see discussions of this in Marcus et al. ( 2011 ) and Levison et al. (j2008 )). 



Proposed models for the formation of the Haumea family have attempted to reproduce 
the family's relatively small velocity dispersion (~ 150 ms~^) and to explain the composi- 
tional and orbital characteristics of the family. However, the orbits of the family members 
have been sculpted by several gigayears of dynamical evolution. In this paper we use nu- 
merical simulations to determine how this orbital evolution affects the dynamical coherence 
of the family. In Section [21 we determine the loss rates for the family, which depend on 
the initial velocity dispersion from the collision, and we determine how the velocity disper- 
sion of the surviving family members is altered over time; from these simulations, we also 
obtain a hard lower limit for the age of the fami l y. In Section [HI we apply these results 



to the family-formation models of iLeinhardt et al\ (120 10|) fa graze-and-mer g e typ e collision 



between two similarly sized, differentiated KBOs) and iSchlichting and Saril (120091 ) (the col- 
lisional disruption of a satellite orbiting Haumea), and we compare the predictions from 
these two formation models to the current observations of the family. Section H] provides a 
summary of our results and conclusions. 



2. Orbital evolution of the Haumea family 

Even though the identified Haumea family members (see Table [1]) have a fairly low ve- 
locity dispersion (Af ~ 150 ms~^), their proper orbital elements span a relatively large range 
in semimajor axis, a, and eccentricity, e, (a range that is typical of classical KBOs), and they 
have atypically large inclinations, z, of ~ 27°. Using the data for their best-fit orbit^, we 
did a 10 Myr numerical simulation to obtain the average values of a, e and i for each family 



^orbit information was taken from the AstDyS website (http://hamilton.dm.unipi.it/astdys 
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member over that time span, and we calculated the correspondin g values of Av (equa- 
tion [p relative to the center-of-mass collision orbit determined by iRaeozzine and Brown 
( 2007); these are listed in T able H] for the family rn e mber s iden tified by Brown et al. ( [20071 ) . 



Schaller and Brownl (|2008[ ) , iRagozzine and Brownl (120071 ) , and ISnodgrass et al\ (120 lOf ) . Be- 



low, we examine the orbital distribution of the known family membe rs to refine the center- 
of-ma ss orbit in light of the additional identified family members since IRagozzine and Brown 
( 120071 ) ■ We use the results of long-term numerical simulations to estimate how much the 
family's orbits have evolved since its formation, and we obtain a hard lower limit on the age 
of the family. 



2.1. Collision center-of-mass orbit and a lower limit on the family's age 

We use the average values of a, e, and i for the nine identified family me mbers (Table [1]) to 



re-calc ulate the center-of-mass collision orbit using the method described by IRagozzine and Brown 



( l2007l ): we minimize the sum of Av for the nine family members while fixing the semimajor 
axis of the center-of-mass orbit at that of Haumea's current orbit, allowing its eccentricity 
and inclination to vary such that Haumea's current Tisserand parameter with respect to 
Neptune (T/v = 2.83) is maintained, and allowing the mean anomaly, M, and the argument 
of pericenter, w, to vary freely; the longitude of ascending node, fi, is ignorable as it does 
not affect the distribution of Av. Figure H] shows the results of this calculation for a range 
of eccentricity and inclination combinations of the collision center-of-mass orbit. The lower 
limit of the shaded region in the figure is the value of the family's average Av found by 
selecting values of the mean anomaly and argument of pericenter that minimize Aw; the 
shaded area shows the range in Av obtained by allowing uj to vary, but still selecting the 
value of M that minimizes Av for each value of u. Parameters along the lower boundary 
of the shaded regions represent collisions occurring very near to the ecliptic plane, while 
parameters along the upper boundary represent collisions at the extreme, off-ecliptic points 
in the orbit (~ 15 — 20 AU above the ecliptic plane). The difference in ave r age A v for 



the different values of is a factor of ~ 2, as noted by IRagozzine and Brownl (120071 ); this 
increase in average Av for off-ecliptic collision points is due to the fact that producing the 
observed family's spread in inclination requires a larger Av at these locations. Because col- 
lisions near the ecliptic are much more probable than off-ecliptic collisions, we choose the 
center-of-mass orbit that minimizes the lower portion of the filled curve in Figure [TJ The 
result is {a,e,i,u,M) = (43.1 AU, 0.12 4, 28.2°, 270°, 76°). Th i s is ve ry similar to the colli- 



sion center-of-mass orbit determined by IRagozzine and Brownl (120071 ): (a, e, i, u, M) = (43.1 



AU, 0.118, 28.2°, 270.8°, 75.7°), indicating that the newer family members do not significantly 
affect the estimate of the collision center-of-mass orbit. The small difference in the eccen- 
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tricity does not much affect the values of Av for the family members (both values of Av are 
listed in Table [1]) because the calculated Av is a fairly flat function of eccentricity within 
± ~ 10% of its minimum. 



In the above calculations, as in iRagozzine and Brown! (120071 ). we assumed a constant 



semimajor axis and conservation of the Tisserand parameter during the chaotic evolution 
of Haumea's orbit. To test the validity of this assumption, we performed numerical simu- 
lations of resonant diffusion within the 12:7 MMR, and we find that Haumea's Tisserand 
parameter can vary by ±0.5%. This is a small variation, but it does affect the allowable 
combinations of e and i for the best-fit center-of-mass orbit. We performed the minimization 
of the sum of Av for the identified family members while allowing e and i to vary inde- 
pendently, and we find a slightly revised best-fit center-of-mass orbit: {a,e,i,u, M) = (43.1 
AU, 0.124, 27.3°, 276°, 70°). This orbit has a Tisserand parameter Tjy = 2.84, which is within 
the range of T/v found in our numerical simulations. If we additionally relax the constraints 
to allow the semimajor axis of the orbit to vary by ±0.15 AU (the approximate range of 
variation in the 12:7 MMR), we find very similar results: {a,e,i,u, M) = (43.1 ± 0.15 
AU, 0.121, 27.3°, 278°, 68°). These alternate minimum Av center-of-mass orbit fits give us 
an estimate of the uncertainties in the orbital parameters: 

(a, e, i, CO, M)cm = (43.1Af/, 0.115 - 0.132, 27 - 28.3°, 270 - 278°, 68 - 76°). 



These orbits all represent collisions near the ecliptic plane, and the IRagozzine and Brown 



(120071 ) estimate falls within the uncertainties. 



We can use the collision center-of-mass orbit to set a lower limit on the age of the 
Haumea family by determining the minimum time necessary for such an orbit to diffuse to 
Haumea's current eccentricity (e = 0.19) in the 12:7 MMR with Neptune. We generated 800 
test particles with initial conditions within the uncertainties of the collision center-of-mass 
orbit found above, randomized their initial mean anomaly, and integrated these for 1 Gyr. 
We find from this simulation that ~ 3% and ~ 6% of the test particles reach Haumea's 
eccentricity by 500 Myr an d 1 Gyr, respectively; this is a slightly lower efficiency than 



Ragozzine and Brownl (120071 ) found from similar simulations (10% had diffused by 1 Gyr), 
but the two results are consistent given that the number of test particles in their simulation 
was only 78. These results allow us to conclude with ~ 95% confidence that the Haumea 
family is older than 1 Gyr. The fastest that any test particle in our simulation diffused to 
Haumea's eccentricity was ~ 100 Myr (Fig. [3]); this indicates that 100 Myr is a strong lower 
limit on the age of the family. 

Another way to estimate the lower limit on the family's age is to examine the precession 
of the orbital planes of the identified family members. The current values of the longitudes 
of ascending node, Q, of the known family members are indistinguishable from a random 
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distribution. Imediately after the family forming collision, the family members will share a 
common line of nodes on the collision center-of-mass orbit plane, but will have different values 
of the other orbital elements. After the collision, the differences in semimajor axes amongst 
the family members cause their orbit planes to precess at slightly different rates. The nodal 
precession rates range from Q4:°Myr~^ to 81°Myr~^ for the family members strongly affected 
by MMRs with Neptune and 69°Myr~^ to 12°Myr-^ for the non-resonant family members 
(as determined from our numerical simulations of their best-fit orbits). Considering the 
differences in these rates, we expect the nodes to be randomized on a 20 Myr timescale for 
the resonant family members and a 100 Myr timescale for the non-resonant family members. 
Both this estimate and the resonant diffusion timescale of Haumea's eccentricity within the 
12:7 MMR indicate that 100 Myr is a hard lower limit on the age of the family. 



2.2. Long-term orbital evolution 

The Haumea family members' range in semimajor axis includes regions affected by 
various MMRs with Neptune. This means that since the time of the collision (at least 100 
Myr ago, but most likely several Gyr ago), the orbital distribution of the family has been 
modified by dynamical evolution, and some family members have probably been removed by 
orbital instabilities. Comparisons of formation models to the current set of observed family 
members must account for this orbital modification and decay of the total population. 

To determine the effects of long-term orbital evolution on the family, we performed a 
suite of eight numerical simulations, each with a cloud of 800 test particles representing family 
members generated with an isotropic di stribution of initial ejec t ion ve locity vectors about the 



collision center-of-mass determined by iRagozzine and Brown! (120071 ): (a, e, i, M) = (42.1 



AU, 0.118, 28.2°, 270.8°, 75.7°). (As discussed in Section [2?T| there is some uncertainty in 
this collision center-of-mass orbit, but the simulation results do not strongly depend on 
the exact values chosen, as all of the allowed range results in test particles spread over 
very similar ranges in a, e, and i.) In each of the eight simulations, we adopted differ- 
ent values of the magnitude, Af , of the initial ejection velocity of the cloud of test par- 
ticles: At; = 50, 100, 150, 400 ms~^, respectively. We then integrated these test parti- 
cle clouds forward in time for 4 Gyr under the influence of the Sun and the four outer 
planets (in their current con figuration), using the symplectic orbital integration method of 



Wisdom and HolmanI (jl99ll ). Any test particles that approached within a Hill sphere of 



Neptune were considered lost, because these are unstable on very short timescales. 

Figure |2] plots the a — e and a — i distributions for two of our simulations; both the initial 
distributions and a snapshot at 3.5 Gyr are shown. In these plots, the test particles are color 
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coded according to their stability (a particle is considered to be unstable if it approaches 
within a Hill sphere of Neptune at any point in the simulation); most of the test particles 
that are unstable are located either near the inner edge of the family (where their initial 
perihelion distances are nearly Neptune crossing) or near the labeled MMRs with Neptune. 
Previous studies have fou nd that this region of the Kuiper belt is strongly affected by even 
fairly high order MMRs Jchiang et alWimi . Ihvkawka and Mukail lioOTI : I Volk and Malhotra 
2OIII ). so it is not surprising that the family is strongly affected as well. The effect of the 
resonances on the unstable test particles is to increase their eccentricity until they become 
Neptune crossing and then have close encounters with Neptune. A few percent of the test 
particles survive in resonance until the end of the si mulation; this result is consistent with 
the resonant fraction found by iLykawka et al\ (120121 ) in a similar study of the evolution of 
the Haumea family. These stable resonant test particles are mostly found in the 3:2 and 
7:4 MMRs; in these cases, the lo ng-lived test particles are additionally stabilized due to 
the Kozai resonance (jKozailll962l ). The Kozai resonance causes a test particle's argument 
of perihelion to librate about 90° which ensures that perihelion occurs well away from the 
echptic plane, protecting the test particle from close encounters with Neptune. 

The fraction of test particles that survive as family members up to 1.5 and 3.5 Gyr 
are listed in Table [2] for each value of Aw; for initial Af of 50 — 200 ms^^ (typical of the 
observed family) 20 — 25% of the family is lost by 3.5 Gyr, which is consistent with the 
Lykawka et al\ (120121 ) results. In addition to the erosion of the population, we also find that 
the dynamical clustering in proper elements grows weaker (and the apparent Af of the test 
particles increases) over the course of the simulations. We determined the Af for each test 
particle by calculating its proper a, e and i (taking the average over the last 50 Myr of the 
simulation) and then allowing the values of the orbital angles to vary until we find the smallest 
difference between the test particle's orbital velocity and the orbital velocity of the collision's 
center-of-mass orbit. We find that the chaotic diffusion in orbital elements induces a spread 
of 50 — 100 ms~^ among the test particles and shifts the average to slightly higher than their 
initial value of Af . Figure|l]shows the distributions of At; at 3.5 Gyr for three simulations in 
which the initial Aw had values of 50 ms~^, 150 ms~^, and 350 ms~^. We conclude that, while 
some individual test particles that are strongly affected by MMRs with Neptune experience 
larger changes in A?;, for the non-resonant Haumea family members, it is likely that the 
value of Af they acquired at the time of the collision was within =b50 — 100 ms~^ of their 
current At;. 
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3. Implications for family formation 



The numerical simulations described above sliow liow long-term orbital evolution will 
affect the Haumea family's dynamical clustering. In this section, we examine the i r nplica tions 
of those results for the pr oposed family formation models of iLeinhardt et al\ (120 lOf ) and 
Schlichting and Saril ( 120091 ) . These models make specific predictions for the mass and velocity 
distribution of the family members immediately following the formation of the family. We 
use our simulations of the family's orbital evolution to evolve the models' predicted mass 
and velocity distributions to the current epoch so we can compare the models' predictions 
to the currently observed family. 



3.1. Overview of proposed formation models 



Brown et all (120071 ) estimated that the Haumea coUisional family was created in a catas- 
trophic collision between a proto-Haumea of radius R ~ 830 km and another KBO with a 
radius of ~ 500 km. However, the low velocity dispersion amongst the observed family mem- 
bers is problematic for such a model. In a catastrophic collision between two large KBOs, 
the velocity dispersion of the family members should be close t o the escape velocity o f the 
largest remnant, Av ~jw,Haumea ~ 900 ms~^ (see discussion in ILeinhardt et al\ (120101 ) and 
Schlichting and Saril ( 120091 ) ): in contrast, the observed family has Aw = 50-300 ms~^. 



To explain the low velocity dispersion of the observed family, ISchlichting and Saril (|2009[ ) 
propose that the family originates from the catastrophic disruption of a satellite orbiting 
Haumea, rather than the disruption of a proto-Haumea. This actually requires two different 
collisions: a collision between a proto-Haumea and another large KBO that creates a large, 
icy satellite and gives Haumea its unusual shape and fast spin, then a subsequent collision 
between the satellite and another KBO. The latter would create a collisional family with 
values of Av close to the escape speed of the satellite rather than of Haumea. Assuming a 
primarily water ice composition, these authors estimate that the disrupted satellite would 
have had a radius R ~ 260 km to account for the mass of Haumea's remaining satellites 
and the rest of the collisional family; the expected At; would be ~ 200 ms~^ for a satellite 
of this size. For a collisional family formed in this way, the authors estimate that the total 
mass of the family at formation would be no more than ~ 5% of the mass of Haumea 
{Mh ~ 4.2 X 10^1 kg). 



Leinhardt et al\ (120101 ) propose an alternative formation mechanism for the family: a 
graze and merge collision event between two similarly sized, radius R ~ 850 km, differentiated 
KBOs. In this scenario, the two impacting bodies merge after the collision, resulting in a 
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very fast rotating object. The family members and satellites are a result of mass shedding 
due to the high spin rate of the merged body (Haumea) rather than being direct impact 
ejecta; this accounts for the low velocity dispersion of the family members. The authors 
ran several collision simulations and found that the resulting family would have a mass of 
~ 4 — 7% Mhi most of which would be found within ~ 300 ms~^ of the collision orbit. 



3.2. Dynamical evolution of the family 

To determine how the dynamical evolution of the family will affect the mass and 
estimates from these formation scenarios, we generate synthetic families for the two models. 
We detail the assumptions for each step below, but the overall procedure is as follows: we first 
sample the size distribution predicted for the coUisional model to generate a list of family 
members and their sizes. We then assign these members an initial Av according to the 
distribution of velocities from the model. This creates a snapshot of the family immediately 
after formation. To account for 3.5 Gyr of orbital evolution, we randomly assign an orbital 
history from one of our eight simulations (Section |2]) to each of our synthetic family members 
(assigning it from the appropriate simulation based on initial At;). Having assigned each 
family member an orbital history, we then take a snapshot of the surviving family members' 
orbital element and mass distributions at 3.5 Gyr. From this we can calculate the expected 
mass vs. At; distribution of the family for each model. 



For the graze-and-merge coUisional model, iLeinhardt et al\ (120101 ) provide plots of the 
cumulative number of coUisional fragments as a function of fragment mass and the cumu- 
lative mass of fragments as a function of At; from each of their high-resolution collision 
simulations (their Figure 3). We use the given total mass of coUisional fragments to convert 
the normalized number of fragments presented in their Figure 3a to an absolute number of 
fragments. To convert the cumulative mass distribution to a cumulative size distribution, 
we assume that all the fragments have a uniform density of 1.15 g cm~^, which is consistent 
with the family being composed of 80% water ice by mass (their Table 2). We construct 
a set of synthetic family members from this size distribution and assign each fragment an 
initial Aw based on the mass vs. velocity distribution (their Figure 3b); each bin in Af is 
filled with randomly selected family members until the specified mass in that bin is reached. 
Based on its initial value of Au, each family member is then randomly assigned a 3.5 Gyr 
orbital history from our numerical simulations; in this way some family members are re- 
moved from the popu lation, and the others diffuse in orbital elements and apparent Ai;. 



Leinhardt et al\ (120 lOf ) detail four different successful graze-and-merge collision simulations 



with slightly varying initial conditions (their simulations 1 through 4). For each of these 
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family forming simulations, we generate 1500 synthetic, evolved families using the procedure 
above. 



Schlichting and Saril ( 120091 ) did not perform collision simulations for the satellite breakup 



model, so we have to make some assumptions about the size and velocity distribution of the 
resulting family. We adopt a differential size distribution for the family 



N{R)dR oc R^-^dR 



(2) 



where R is the radius of the collision fragments; we adopt values of P in the range 4.5 



5.5, consistent with typical catastrophic collision simulations ( iLeinhardt and StewartI 12012 



Marcus et a/.ll201ll ). Assuming a total family mass of 0.04 — 0.05Mj:/, we generate a synthetic 
family from this size distribution, with fragments in the size range 50 km < R < 150 km (the 
same size range as the graze-and-merge simulations, for ease of comparison). We use the 
same density for the family members as for the graze-and-merge formation scenario to convert 
radius to mass. For the ejection velocity of the fragments, we adopt a normal distribution 
with mean 200 ms~^ and standard deviation 50 ms~^. We generate 1500 synthetic families 
using these assumptions and dynamically evolve the family members for 3.5 Gyr in the same 
way as described for the graze-and-merge model. 



Figures [5] and [6] show our results for the ILeinhardt et all (120101 ) and lSchlichting and Sari 



(120091 ) formation models, respectively: we plot the average cumulative mass of the synthetic 
families as a function of their apparent Av at t = 3.5 Gyr; the la uncertainties are also 
indicated. Our calculations find that in the graze-and-merge model, the current evolved 
Haumea family should have a total mass of 0.045 ± 0.01 M^, of which ~ 0.02 Mh should 
be found at Av < 150 ms~^. The currently o bserved farnily (^in cluding Haumea's satellites) 
is estimated to have a mass of ~ 0.017 Mh (ICook et a/.ll201ll ). which accounts for ~ 85% 
of the mass tha t is ex pected to be found within 150 ms~^ of the collision center for the 
Leinhardt et all (120101 ) formation model. The model predicts that there should be twice as 
much mass (0.035 ± 0.01 Mh) in family members at larger velocities. The satellite breakup 
model predicts a surviving family of ~ 0.035 Mh, mostly in the Av = 100-300 ms~^ range, 
indicating that the known family members account for ~ 50% of the expected mass of the 
family. 



3.3. Comparison with observations 

The known Haumea family members are not an observationally complete population; 
to compare the evolved synthetic families from Section 13.21 to the observed family, we must 
estimate how many of our synthetic family members would be within the observable ap- 
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parent magnitude and ecliptic latitude range of observational surveys that have detected 
the Haumea family members. For each of our synthetic families (obtained in Section |3T2|) . 
we take a snapshot of the instantaneous orbital elements of the family members at t = 3.5 
Gyr from their assigned orbital histories. This snapshot allows us to calculate the helio- 
centric distance and ecliptic latitude for each synthetic family member. To calculate the 
apparent magnitude, we use the heliocentric distance and the assigned size (as described in 
Section [3l2|) . but we also need to r nake some assumptions ab out t he albedos of the f amily 
members. Haumea's albedo is 0.8 ( iLacerda and JewittI 120071 ). and lEUiot et al\ (120101 ) have 
measured an albedo of 0.88 for the next brightest family member (2003 TX300), but albedos 
have not been measured for the other family members. Based on the light curves obtained for 
five of the know n family members and th e light curves of other icy solar system bodies with 
known albedos, iRabinowitz et al\ (120081 ) argue that the Haumea family members' albedos 
are likely in the range 0.3 — 1.4. For our synthetic family, we adopt an albedo distribution 
with an average of 0.8 (Haumea's albedo) and a uniform spread of ±0.2. (We discuss below 
how this assumption might affect the results of our comparison.) Given this albedo assump- 
tion, we calculate the apparent magnitudes and ecliptic latitudes for each synthetic family, 
and we use the resulting distributions to calculate the number of objects as a function of Av 
for each synthetic family that would be detected by an observational survey. 

For the observati onal comparison , we u se the results of the Palomar distant solar system 
survey conducted by ISchwamb et all (120101 ). This was a wide-field survey of ~ 12000 deg^ 
down to a limiting magnitude of m,, ^ 21.3, detecting 52 KBOs and Centaurs (27 previ- 
ously known objects, and 25 new ones), including 4 of the previously identified Haumea 
family member s . Th e presence of so many known KBOs in their survey fields allowed 
Schwamb et al\ ( I2OIOI ) to estimate that their detection efficiency was ~ 65% down to nir — 



21.3 for the known population (see their Figure 3). They also provide a plot of the survey's 
fractional sky coverage as a function of ecliptic latitude (see their Figure 4); they covered 
approximately 50% of the sky ±30° from the ecliptic. From this information, we can esti- 
mate the detection probability for an object in our synthetic families based on its apparent 
magnitude and ecliptic latitude. We use this detection probability to determine how many 
of our synthetic coUisional family members could have been detected in this survey. 



An important note here is that the ISchwamb et all ( 120101 ) survey was not capable of 
spectrally identifying family members. They can only say that there are four previously 
identified Haumea family members within their detections (listed in their Table 2), but it is 
possible that additional, unidentified Haumea family members are present within their survey 
detections. To examine this possibility, we calculated Av for each of their listed detections 
and found two additional objects within 500 ms~^ of the collisio n center-of-mass orbi t . One 
of these objects, 2004 SBgo {Av ^ 350 ms~^), was observed by ISchaller and Brownl (|2008[ ) 
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and found not to have the water ice spectral feature characteristic of all the other identified 
Haumea family members (its surface spectrum is consistent with no water ice being present). 
The other object, 2008 AP129, has a Ai; f» 140 ms~^, suggestive of a dynamical association 
with the Haumea family, bu t the object shows only a moderate amount of water ice absorp- 



tion in its surface spectrum (IBrown et a/.ll2012l ). with a water ice fraction substantially lower 



than measured for the known family members. If we allow for the possibility that 2008 AP129 



i s a w ater ice poor family member, this means that it is possible that the ISchwamb et al. 



(120101 ) survey detected as many as five Haumea family members, but this number is not 



likely to be larger. 

Figures [7] and [8] show the number of synthetic family members as a fu nction of Af for 



each of the two formation scenarios that would have been detected by the ISchwamb et al 



(120101 ) survey; also shown for reference are the Haumea family members actually detected 



in that survey. 

For the satellite breakup model (Figure [8]), the actual number of family members de- 
tected falls within the la range of total detections predicted by the model, but the observed 
values of Av are lower than the predicted values. Almost none of the synthetic families 
match the observations by producing 4 or 5 detections all with Av < 150 ms~^. The values 
of Av for some of the observed family members could be larger than the values in Table [1] 
though, because of the uncertainty in the orbital angles of the collision orbit. To constrain 
the collision orbit and calculate the minimum Av for the known family members, we assumed 
that the collision took place near the ecliptic plane (where collision probabilities are highest); 
the family mem bers' ejection velocities from the collision could be larger than this minimum 



value, although iRagozzine and Brownl ( l2007l ) argue that the correction factor is likely to be 
~ 2 or lower unless the ejection of fragments was highly anisotropic (see also our discussion 
of this in Section |2?T1) . If we allow for a factor of two correction to the Aw estimates for the 
real family members, we can extend the Av of the survey's detected family members out 
to ~ 250 ms~^. Using this increased allowable range of Av, 18% of the synthetic families 
satisfactorily reproduce the observations, making the satellite breakup model statistically 
consistent with the observations. 

The graze-and-merge model predicts that, on average, the survey should have detected 
8 family members. This is larger than the 4 or 5 real detections, but the actual detections 
fall within the la uncertainties of the synthetic families. Just as in the satellite breakup 
model, the real detections fall at significantly lower Av than predicted by the graze-and- 
merge model. Very few of the synthetic families result in all the detectable family members 
falling below 150 ms~^, and all of these cases result in too few detections to be consistent 
with the observations. If we increase the allowed maximum Av to 250 ms~^ and require 



- 13 - 



the synthetic famihes to match the number of real detections (4 or 5), we find that only 
4% of the synthetic families reproduce the observations, indicating that the observations 
are not a typical outcome of the graze-and-merge model. We did, however, make a number 
of assumptions when generating our synthetic families, so we examine how these might be 
altered to bring the model into closer agreement with the observations. 



One assumption we made in the creation of the families in Section 13.21 was that there 
was no relationship between a fragment's mass and its Av from the collision center. Our only 



constraint was that the binned mass vs. Av matched the outcome of the Leinhardt et al. 



( 120101 ) simulations. If there is a correlation between fragment mass and Av such that the 
higher Aw family members are, on average, smaller than the lower Av members, this might 
account for the lack of observed family members at large Av. To test this we impose a 
relationship between fragment mass, m, and Av such that, averaged over all the family 
members, 



Avirn] 



oc m 



(3) 



while the total mass o f all fragments within a given range of Av is still constrained by the 

Leinhardt et al\ (120101 ) simulation results. This power law relationship has been seen in some 

laboratory impact experiments, with the expone nt typically being /c < 1/6 (INakamura and Fujiwara 
199ll : iHolsapple et a/.ll2002l : iGiblin et a/.l 120041 ). Taking /c = 1/6 and generating a new set 
of synthetic families for the graze-and-merge collision, the percentage of synthetic families 
with 4 or 5 detections all with Aw < 250 ms~^ increases to 8%. The percentage of synthetic 
families with 4 or 5 detections and Av < 150 ms~^ is still negligible; even if we increase the 
value of k to 1/4 (which is not a likely value), we still see < 1% of the synthetic families 
resulting in 4 or 5 detections below 150 ms~^. 

Another easily altered assumption is that of the albedos for the family members. We 
assumed albedos in the range 0.6 — 1.0, but if we assume systematically lower albedos, in the 
range 0.4 — 0.6 (still consistent with their icy composition), we decrease the average predicted 
number of detections for the graze-and-merge model from ~ 8 to ~ 4, also decreasing the 
number of predicted detections at large Av. With this albedo assumption, nearly half of the 
synthetic families result in no detections at Av > 250 ms~^, and the percentage of families 
with 4 or 5 detections and Av < 250 ms~^ is 11%. Lowering the albedos further does 
not increase the percentage of families that agree with the observations. If we assume the 
above albedo range in combination with the k = 1/6 mass- Aw relationship, the percentage 
of matching families is ~ 12%. Given that at least one family member besides Haumea 
has been shown to have a very bright albedo ( lElliot et a/.l 120101 ). assuming systematically 
lower albedos for most of the family members might not be the most likely solution to the 
discrepancy between the observations and the models. However, it cannot be ruled out until 
the albedos of some of the smaller family members have been measured. 
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Depending on the assumptions about albedos and mass-velocity relationships, and al- 
lowing for the uncertainty in the known family's Av, we find that 4 — 12% of the synthetic 
families generat ed by the graze-and-me rge collision model reproduce the Haumea family, as 
observed by the ISchwamb et al\ (120 10[ ) survey. The actual observed family is at systemati- 
cally lower Av than the model predicts, but given the relatively small set of observed family 
members, the model is still consistent the observations. 



4. Discussion and Conclusions 

After accounting for 3.5 Gyr of dynamical evolution and accounting for the observational 
incompleteness of the k nown family, both of the proposed Haumea fami ly formation models 



we have examined here (ILeinhardt et a/.ll2010t ISchlichting and Sarill2009l ) are consistent with 



the total number of observed family members. There is, however, a significant difference 
between the observed values of Av and those predicted in the formation models; almost 
none 1%) of the synthetic families we generated for either formation scenario account for 
the observations of family members that all fall within 150 ms~^ of the collision center. The 
only way we find to make the velocity distributions for the formation scenarios consistent with 
the observations is to allow that the actual ejection velocities of some of the known family 
members are a factor of ~ 2 larger than the calculated minimum values. This adjustment of 
the observed Av falls within the uncertainties for their calculation of the collision center-of- 
mass orbit (see our discussion in Section [211!) . But even with this adjustment, only 10—20% of 
the synthetic families reproduce the observed family. This is a reasonable level of agreement 
given the uncertainties in the collision models and the small number of observed family 
members, but it is still interesting that so few family members have been identified at large 
Av. In section [23] we explored how the assumptions we made about the albedos of the family 
members and possible correlations between fragment mass and ejection velocity could change 
the predictions of the graze-and-merge collision model. We find that it is difficult to alter 
these assumptions enough to obtain a better than 10-15% agreement with the observations. 
For either formation scenario, we find that there should be an additional ~ 0.01 — 0.03 Mh 
of family members that have yet to be identified, i.e., at least as many as those already 
detected. If, as we discover these additional family members, the distribution of Av remains 
too heavily weighted toward the low-end, the formation models will have to be reconsidered. 

Recent spectroscopic and photometric surveys of the Kuiper belt desi gned to detect wa- 



ter ice have not identified any large Aw, ice- rich Haumea family members. iFraser and Brown 



(120121 ) performed a photometric study of 120 objects from the major dynamical classes of the 



Kuiper belt using the Hubble Space Telescope (HST) and did not find any additional Haumea 
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family members. iBenecchi et all (120111 ) also performed HST photometry of a large sample 
of Kuiper belt objects and failed to identify any new ice-rich family members. Ground based 
studies searc hing for water ice in the Kuiper belt have also failed to detect hi gher Av fam- 



i ly me mbers; iBrown et all ( 120121 ) detected no additional family members, and iTrujillo et al 



( I2OIII ) found only one new member, located near the dynamical core of the family. The 
extent of these surveys suggests that if there were ice-rich Haumea family members spread 
at large Aw throughout the Kuiper belt, we likely should have already identified some of 
them. This is consistent with our comparisons of the proposed collisional models to the 
observed Haumea family; figure [9] shows the distribution of simulated detections for a subset 
of our dynamically evolved graze-and-merge collisional families (see Section 13. 3p compared 
to the known Haumea family. The simulated detections span a larger parameter range of 
the K uiper belt (due to the presence of higher Aw family members) than the actual detec- 
tions; iLykawka et al\ ( 120121 ) also found that simulated Haumea families with Aw consistent 
with existing collisional formation models would occupy a larger portion of the Kuiper belt 
than currently observed. In contrast to the expected Au distributions from the collisional 
models, we find that, accounting for the time evolution of the a-e-i and Aw distributions in 
our simulations, all of the known family members are consistent with initial Af < 100 ms~^. 



One possible way to modify the collision models was suggested by ICook et al\ ( 1201 ll ). 
who argue that perhaps the proto-Haumea was only partially differentiated. This would allow 
for some of the collisional fragments to be rocky rather than primarily water ice, which would 
mean that they might not show the water ice spectral feature that has thus far been used 
as the only secure way to identify a family member; perhaps the dynamically nearby KBO 
2008 AP129, which has a Af of only 140 ms~^ but shows less water ice absorption than the 
accepted family members (see Section 13. 3p represents a new class of rockier family members. 
A rockier composition could also make the fragments darker and therefore more difficult 
to detect at all. It has been similarly suggested that surface inhomogeneities on the proto- 
Haumea could result in collisional fra gments with different comp ositional characteristics from 
those of the known family members ( ISchaller and Brownl l2008l ) . It is unclear why, in either 
of these situations, there would be a preference for the less icy fragments to be dispersed 
at large At;, but if the composition of the target and impactor are substantially different 
from those assumed in the collision simulations, that could affect the entire Av distribution. 
Another possibility is that the formation models have not adequately accounted for collisional 
evolution amongst the ejected fragments themselves, and that this could alter the family's 
siz e and/or velocity dist ribution in a significant way. Collision simulations like the ones 



m 



Leinhardt et al\ ( I2OIOI ) are computationally very expensive and they follow the family's 



evolution for only a few thousand spin periods of the primary (a few hundred days in total), 
so the model's final size and velocity distributions are not fully evolved. 
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The presence or absence of higher Av family members with future additions to the 
set of observed Haumea family members will determine if any of these modified scenarios 
should be considered, or if there is another c ollision al model that could better explain the 
family. For either the iLeinhardt et al. J2010h or the ISchlichting and Saril (120091 ) models to 
be consistent with observations, we should find several higher Av family members among 
any new identifications. 

In summary, our study of the long-term dynamical evolution of the Haumea family leads 
to the followings conclusions. 



2. 



4. 



The Haumea family is at least 100 Myr old. This estimate is based on the timescale 
to randomize the nodal longitudes of the orbital planes of the family members, as well 
the timescale for chaotic evolution of Haumea's eccentricity in the 12:7 MMR with 
Neptune. From the chaotic diffusion of Haumea's eccentricity, we can conclude with 
95% confidence that the family is older than 1 Gyr. 

For initial ejection velocities. Aw, in the range 50 — 400 ms~^, 20 — 45% of original 
Haumea family members are lost due to close encounters with Neptune over 3.5 Gyr. 
Most of this loss occurs at the inner edge of the family (interior to ~ 41 AU) and near 
the locations of MMRs with Neptune. A few percent of the surviving Haumea family 
members are expected to be found in MMRs with Neptune. The 3:2 and 7:4 MMRs 
are the most likely of the resonances to contain surviving members. 

Within the population of surviving and potentially recognizable family members, chaotic 
diffusion in orbital elements over 3.5 Gyr introduces a 50 — 100 ms~^ spread in the 
apparent velocities of the family relative to the collision center-of-mass orbit, with the 
average Av increasing slightly over time. 

A ccounting for long-terra dynamical evolution to the graze-and-merge collision model 
of lLeinhardt et al\ (120101 ). we find that the currently observed family represents > 85% 
of the expected family mass within 150 ms~^ of the collision center, but an additional 
0.035 ± 0.01 Mh (about twice the mass of the known family) remains to be identified 
at larger Av. Accounting for observational incompleteness, the ILeinhardt et al\ (120 lOf ) 
model is consistent with the observations at the ~ 10% confidence level. 



For the satellite breakup model of ISchlichting and Saril ( 120091 ) . we find that the cur- 
rently observed family accounts for ~ 50% of the expected mass of the family. Most of 
the remaining mass should be found at Av > 150 ms~^. Accounting for observational 
incompleteness, the satellite breakup model is consistent with the observations at the 
~ 20% confidence level. 
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6. Both formation models predict more family members at large Av than are currently 
observed (even allowing for a factor of ~ 2 higher values of Av for the known family 
members due to the uncertainty in estimates of the collision center-of-mass orbit). 
If additional Haumea family members are identified and continue to have low Av 
(< 200 ms~^), new formation models (or modifications to the existing models) will 
have to be considered. 



This research was supported by grant no. NNX08AQ65G from NASA's Outer Planets 
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Fig. 1. — Average Av for the nine identified Haumea family members as a function of the 
eccentricity and inchnation of the coUision center-of-mass orbit. The shaded regions show 
the range of Av over all values of the argument of pericenter (u) for the center-of-mass orbit 
(assuming the value of the mean anomaly that minimizes Av for each value of u). 



Table 1. Known Haumea Family Members 



Designation 


a (AU) 


e 


i (deg) 


Avi (ms ^) 


Av2 (nis^^)-^ 


absolute magnitude 


Haumea (2003 ELgi) 


43.10 


0.19 


26.8 


320 


342 


0.27 


1995 SM55 


41.85 


0.10 


26.7 


158 


174 


6.3 


1996 TOee 


43.34 


0.12 


28.0 


25 


46 


4.5 


1999 OY3 


43.92 


0.17 


25.8 


293 


272 


6.76 


2002 TX300 


43.30 


0.13 


27.0 


107 


97 


3.09 


2003 OP32 


43.25 


0.10 


27.0 


120 


138 


4.2 


2003 UZ117 


44.28 


0.13 


28.0 


67 


60 


5.2 


2005 CB79 


43.40 


0.13 


27.2 


111 


94 


5.4 


2005 RR43 


43.28 


0.13 


27.1 


113 


98 


4.0 


2003 SQ317 


42.68 


0.09 


28.1 


144 


172 


6.3 



Note. — Avi is calculated using the iRagozzine and Brownl (|2007( ) center-of-mass orbit and Av2 is 
calculated using our updated center-of-mass orbit (see section [^TT|) . When calculating Av, a, e, and i for 
the family members are fixed at their current values, but the orbital angles are allowed to vary freely. 
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Fig. 2. — Proper eccentricity and semimajor axes for the observed family members (green circles) and for 
sets of test particles with an isotropic Aw = 150 ms^^ (top four panels) and Aw = 350 ms~^ (bottom four 
panels) from the collision ccntcr-of-mass orbit (red triangle). Black points indicate test particles that are 
stable over 3.5 Gyr; blue points indicate test particles that are unstable and are removed via close encounters 
with Neptune within 3.5 Gyr. The locations of various MMRs with Neptune are indicated. Note that the 
scaling of the x- and y-axes is different in the top four and bottom four panels. The two values of Aw 
are shown as examples of the evolution of the family. The simulations for Aw = 50, 100, 200, 250, 300 and 
400 ms^^ are not shown. 
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Fig. 3. — Eccentricity time evolution for 800 clones of the best-fit collision center-of-mass 
orbit. The grey dots are for particles which have not evolved to Haumea's current eccentricity 
by the end of the 1 Gyr simulation (94% of the cloned orbits). The black dots are the orbits 
which do reach Haumea's eccentricity (6% of the cloned orbits). The minimum time required 
to diffuse to Haumea's eccentricity is just over 100 Myr. 
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Fig. 4. — The distribution of relative velocity, Af , for the surviving test particles at 3.5 
Gyr in the two numerical simulations shown in Figure |2] as well as the simulation with 
initial Av = 50 ms~^; in each simulation we started with a cloud of particles having the 
same magnitude of Av (of 50, 150, and 350 ms~^), but isotropically distributed in direction, 
relative to the velocity of the collision center-of-mass orbit. 
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Fig. 5. — Average cumulative mass (in units of Haumea's mass) as a function of apparent 
Av aX t = 3.5Gyr (solid line) with la variations (dashed lines) for 6000 synthetic graze- 
and- merge families (see section [3l2|) . For comparison, the total mass of the known Haumea 
family is estimated to be ~ 0.017 Mh and is located mostly with Av < 150 ms~^. 
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Fig. 6. — Average cumulative mass (in units of Haumea's mass) as a function of apparent 
Av at t = 3.5Gyr (solid line) with la variations (dashed lines) for 1500 synthetic satellite- 
breakup families families (see section 13. 2p . 
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Fig. 7. — Cumul ative velocity d i stribu tion of synthetic graze-and-merge family members 
observable by the ISchwamb et al\ (120101 ) distant solar system survey compared to the distri- 
bution of Haumea family members actually detected (see section [3l3!) . 



<i 

A 



10 



observed 

average of syntlietic families 
1 a variation 




400 



500 



A V (ms" ) 



Fig. 8. — The same as for Figure [7] but for the satellite-breakup formation model. 
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Fig. 9. — Eccentricity vs. semimajor axis (left panel) and inclination vs. semimajor axis 
(right panel) for the simulated detections (black dots) from a representative subset of the 
synthetic graze-and-merge collisional families. The observed family is shown in grey. 



Table 2. Fraction of simulated test particles that survive to 1.5 and 3.5 Gyr 



Av (ms "'") 


t = 1.5 Gyr 


t = 3.5 Gyr 


50 


0.89 


0.81 


100 


0.87 


0.77 


150 


0.88 


0.80 


200 


0.79 


0.74 


250 


0.75 


0.67 


300 


0.70 


0.64 


350 


0.65 


0.61 


400 


0.63 


0.57 



